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Abstract. 

The RHESSI experiment uses rotational modulation for x- and gamma ray imaging 
of solar eruptions. In order to disentangle rotational modulation from intrinsic time 
variation, an unbiased linear estimator for the spatially integrated photon flux is 
proposed. The estimator mimics a flat instrumental response under a gaussian prior, 
with achievable flatness depending on the counting noise. The amount of regularization 
is primarily given by the modulation-to-Poisson levels of fluctuations, and is only 
weakly affected by the Bayesian prior. Monte Carlo simulations demonstrate that 
the mean relative error of the estimator reaches the Poisson limit, and real-data 
applications are shown. 

PACS numbers: 95.55.Ev, 95.55.Ka, 95.75.Wx, 95.75.Pq 
1. Introduction 

At present, imaging of hard x-rays (HXR) and Gamma rays (GR) is only feasible by 
selective absorption using different kinds of masks. An elegant and economic variant 
uses rotation collimators (Schnopper 1968, Willmore 1970, Skinner and Ponman 1995) 
where a pair of absorbing grids rotates between the true scene and a spatially non- 
resolving detector (Fig. [T]top). Depending on whether a source is behind or between 
the grid bars, the observed flux is low or large. As rotation progresses, the true scene 
becomes thus encoded in a temporal modulation of the observed HXR/GR flux (Fig. 
n bottom). In this process, the modulation frequency is not constant but varies with 
time: glancing passages of the grid bars produce slow modulation, and rippling passages 
yield fast modulation. The absolute value of the modulation frequency also depends on 
the offset of the source from the rotation axis, and the modulation amplitude depends 
on the source size compared to the grid period, in such a way that the amplitude is 
largest for a point source, and tends to zero if the source size exceeds the grid period. 
Altogether, this results in a characteristic time series, from which the true scene can 
be reconstructed by suitable inverse methods (Skinner 1979, Prince et al 1988, Skinner 
and Ponman 1995, Hurford et al 2002a). While such methods usually assume that the 
true scene does not depend on time, the present article deals with the complementary 
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problem of estimating the true time dependence of the (spatially integrated) scene, 
and distinguishing it from rotation modulation. This is called here the ^demodulation'' 
problem. 

The data which we envisage are obtained by the Reuven Ramaty High-Energy Solar 
Spectroscopic Imager (RHESSI, Lin et al 2002; Zehnder et al 2003). This solar-dedicated 
space mission uses rotation modulation for HXR/GR imaging of solar eruptions. The 
RHESSI instrument has 9 pairs of aligned HXR/GR-absorbing grids (9 'subcollimators') 
which are fixed on the rotating spacecraft. Different subcollimators have different 
grid spacings, which produce sinusoidal spatial transmission patterns (see Section El 
for details). Behind each subcoUimator there is a germanium detector which records 
the arrival time and energy (3 keV - 17 MeV) of the incoming HXR/GR photons. From 
a statistical point of view, the photon arrival times form a Poisson process with non- 
constant intensity. The change in intensity is due to rotation modulation, but may also 
have contributions from the time dependence of the true scene. 

The latter is highly interesting from a solar physics point of view, because it is 
related to particle acceleration in impulsive solar eruptions. The existence of temporal 
fluctuations of the true HXR scene down to time scales of hundred milliseconds has been 
confirmed by earlier observations (Dennis 1985, Machado et al 1993, Aschwanden et al 
1993). However, these time scales interfere with the RHESSI rotation modulation which 
occurs on time scales of milliseconds to seconds. In order to disentangle time dependence 
of the true scene from rotation modulation the RHESSI data must be demodulated. The 
most general outcome of this process would be the true scene as a function of both space 
and time. But the information contained in the photon counts is rather limited, and 
must be shared between spatial and temporal degrees of freedom of the true scene. 
We therefore restrict ourselves here to the simpler problem of estimating the spatially 
integrated true scene as a function of time. The motivation for this stems also from 
the wish to compare RHESSI HXR data with spatially integrated broadband radio 
observations. Both types of radiation are emitted by electrons of comparable energy, 
and there is some controversy in the field whether or not the two populations actually 
agree. The actual radio observations have a time resolution of 100 milliseconds, and 
collect radiation from the whole sun. 

The demodulation of RHESSI data is, in general, an inverse problem. Temporal 
and spatial variations of the true scene are entangled by the observation method. An 
exception only arises if the time scale of interest exceeds the RHESSI rotation period Ts 
~ 4s, so that demodulation can be replaced by a running average over Ts- Otherwise, 
a good estimator for the spatially integrated scene requires some a priori information 
on its spatial structure in order to outweigh or damp the rotational modulation. Such 
information may either come from independent observations, or from RHESSI itself. In 
the latter case, one may use standard RHESSI imaging techniques (Hurford et al 2002a) 
to obtain an estimate of the time-averaged true scene. These techniques assume that 
the true scene is independent of time. Under this assumption, the true scene can be 
estimated after half a spacecraft rotation, when all possible grid orientations are cycled. 
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and higher- quality results arise from multiples of T5/2. 

The aim of the present paper is to resolve time scales <^ Ts/2. We do not attempt 
here a fully general solution where the spatially true integrated scene is found at each 
point of time. Instead, we will assume that the true scene can be considered as piecewise 
constant during time intervals r of order of 100 ms. Such intervals are shorter than 
most previously reported time scales, and they also agree with the time resolution of 
the comparative radio observations mentioned above. The present article describes an 
efficient linear estimator for the spatially- and r-integrated scene. The construction of 
the estimator follows the lines of classical Wiener filtering (Rybicki and Press 1992), 
but operates in time (not frequency) space. 

The paper is organized as follows. Section|21summarizes some relevant technicalities 
of the RHESSI instrument. Section IHl presents our inverse theory approach including 
counting statistics (Sect. IH.lj) and a priori information (Sect. 13. 2j) . The resulting 
estimator is then discussed in Sections 14.11 - and its performance (Sect. 14. 4p and 
robustness against violation of the prior assumptions (Sect. 14. 5j) are explored by Monte 
Carlo simulations. Section 14.61 shows the algorithm at work on real data from a solar 
eruption. See Tab. Qfor an overview of notation. 
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= w^c^, estimator for b [ct] 



Table 1. Notation. Angular position x = {x,y) is measured in locally Cartesian 
heliocentric coordinates, and the dependence on x is therefore also termed 'spatial'. 
One arc second (1") corresponds to 700 km on the solar disc. 

2. The RHESSI instrument 

We start with a brief description of the RHESSI response to HXRs and GRs, which 
defines the 'forward' problem of converting the true scene into the observed counts. We 
shall only consider photons out of a fixed energy band, average all energy-dependent 
quantities over that energy band, and omit the energy dependence in the notation. 
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2.1. Modulation patterns 

The instantaneous transmission probabilities of RHESSI's subcoUimators, as a function 
of photon incidence direction, are called the modulation patterns. They may be 
visualized as the grids' shadow on the sky plane if the detectors were operated in 
a transmitting mode. When expressed in terms of heliocentric cartesian coordinates 
X = [x, y) (over the limited RHESSI field of view, cartesian and angular coordinates are 
equivalent), the modulation patterns can be approximated by 

Mi(x, t) = aoi(t) + aii{t) cos (ki{t) ■ (x - P(t)) + ^^(t)) . (1) 

Above, the kj(t), i = 1..9, are the grid vectors which rotate clockwise with the 
spacecraft rotation period T5. All grid periods (27r/|A;*| = 2.6 ■ 3*^^ asec) are small 
compared to the solar diameter (1920") and to the fields of view of the individual 
subcoUimators (3600" ... 27.000"). The coefficients (a„i(t), \E'j(t)) describe the internal 
shadowing of the grids; they depend on photon energy and only weakly on time as long 
as the source is in the central part of the field of view, which is the case for solar sources. 
The vector P{t) is the imaging (optical) axis. It generally varies with time, since it is 
not aligned with the rotation axis, which, in turn, may differ from the inertial axis. As 
a consequence, P(t) traces out a relatively complicated orbit on the solar disc, which is 
continuously monitored by the onboard aspect systems (Fivian et al 2002, Hurford et 
al 2002b), and the details of which need not concern us here. All vectors x, ]<.{t) and 
P(t) ly in the solar plane. 

2.2. Onboard data reduction 

In order to avoid detector and telemetry saturation, mechanical attenuators can be 
inserted into the optical path, and photon counts can be decimated by a clocked veto: 
during one binary microsecond (Ib/is = 2~^°s), all events are accepted, during the next 
{Nfi — l)b/is, they are all rejected. Both reduction methods preserve the statistical 
independence of the photon arriving times, and can be absorbed in the definition of the 
modulation patterns. 

2.3. Observational artifacts 

While Equation represents the ideal instrumental response, the real observations 
suffer from several non-idealities such as detector saturation at high count rates, and 
sporadic breakdown of the whole detection chain. The latter is presumably caused by 
cosmic rays (Smith et al 2002); the resulting data gaps occur at random about once 
per second, and have durations of milliseconds to seconds (see Fig. Q bottom for an 
example). All non-idealities are combined on-ground in livetime measures < Li{t) < 
1, which estimate the operational fraction of each detector in a given time bin. The 
effect of livetime is taken into account by multiplying the modulation patterns Mj(x, t) 
by Li(t). Times with L < 0.3 are discarded, which prevents detector saturation and 
redundant or undefined {L = 0) numerical operations. 
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3. Inverse Method 

3.1. Principle and treatment of measurement errors 

Let -B(x, t) denote the solar brightness distribution at position x and time t, and let the 
RHESSI counts be grouped in time bins G r where = (z, j) labels subcoUimator and 
time. The set {//} may comprise all or only some of the 9 subcoUimators, and different 
sub collimators may have different time bins. The observed counts in different bins 
are supposed to be statistically independent Poisson variates with expectation values 

[ dt [ dxMi(^^)ix,t)B{x,t), (2) 
Ja^ J 

where we have written to extract the subcoUimator index out of fi. Non-solar 
background is neglected. The goal is to estimate spatially integrated true scene 

dt j dxB{x,t) (3) 

in the time interval r, as would be observed in the absence of the subcoUimators 
(M((x, t) = 1). The estimator b for b is searched in the linear form 

i) = Wf^c^ (summation convention). (4) 

In the sequel we shall choose to make b efficient, i.e., unbiased and minimum- 
variance among all unbiased estimators. To this end, and for the sake of a concise 
notation, we first introduce the expectation operator £^ of a general function X{b, c) 
of the true b and of the observed counts c = {c^}. This is done within a Bayesian 
framework where each true scene B has assigned a prior probability P{B), so that 
P{c,B) = P{c\B)P{B). One thus has 

SX{b, c) = / dP{B) n ^(^' ^) (5) 

c=o 

^ ? ' 

where and b depend on B according to Equations (j21) and 0. The expectation 
operator ^ involves an average over the counting statistics, followed by an average over 
the prior scenes, and we write £ = £b £c\b to stress this sequence. If a quantity X is 
independent of the observed counts c then £ reduces to £b- 

So far, £-Q is still a formal object. Let us accept this for the moment, and proceed 
with the design of the weights w = {w^}. A first natural request is that the estimator 
(jH is unbiased, £{b) = £{b), requiring that 

w^£B{X^)=£B{b). (6) 

This condition locates w along the direction of £b{X). In order to find a 
unique solution w we additionally require the variance £{b — b)"^ to be minimal, 
and verify a posteriori that this condition is sufficient to ensure uniqueness of the 
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solution. Unbiasedness and minimum-variance constraints are combined by minimizing 
S(b — b)'^+(£ (b — b) with respect to w, where the Lagrange multipher ( must be adjusted 
to fulfill Equation (jHl). The resulting minimum condition is 

S{Cf,c,y)wu - 8{bcu) + C^{cu) = , (7) 

or, after performing the average £c\b over counting statistics, 

Sb{KK + K5^u)w, - SB{bX,) + CSb{Xu) = . (8) 

Here it was used that Sc\b{c^Cu) = Xfj.Xu + X^6^i, for independent counts c^. One 
may now see that Equation (jH)) has indeed a unique solution w if only > 0. 
Under this condition, the matrix A^^, = X^X^ + X^6^u is strictly positive definite: 
(x. Ax) = (x • A)^ + X^xfi^ > for all x 7^ 0. Since the prior average £b represents 
a sum of positively weighted A's, the matrix £bA is also positive definite and therefore 
invertible - indeed, the condition number of can not exceed the condition number 
of A t 

Let us add at this point an interpretation of the matrix A. For any fixed realization 
of the true scene B, the fiuctuations of the observed counts have two causes: 
instrumental modulation and counting noise. These two causes give rise to the two 
contributions A^Aj^ and A^5^j^ in the matrix A^j,^. At large count rates (A^ ^ 1), A^j, 
is dominated by A^A,y, and at small count rates (A^i < 1) it is dominated by X^6^i^. 
The condition number of A^,^ increases with increasing count rates and is bounded from 
above by (| Ap+max(A^))/ min(A^)|. Thus, higher count rates allow weaker-conditioned 
^ij.u- We may interpret this by saying that the counting noise regularizes the problem 
of inverting modulation, and that the amount of regularization is given by the level of 
Poisson fiuctuations compared to the modulation amplitude. 

3.2. Choice of the prior 

We turn now to the definition of Sb (Eq. El). Formally, the true brightness distribution 
B{'x,t) is an element of a function space, and its a priori probability measure 
dP{B) is a functional. In order to avoid technical complications while retaining the 
basic probabilistic features, we make the following simplifying assumptions: the true 
brightness distribution 5(x, t) is concentrated in a (known) spatial prior, where it may 
have (unknown) substructures which do not vary significantly during the time interval 
r for which the unmodulated counts are to be estimated. Thus i?(x, t) ^ -B(x) for 
t E T. Now we recall that -B(x) enters the problem only upon weighting with the 
modulation paterns (Eq. E)). Therefore the detailed structure of -B(x) on scales which 
are small compared to the period of the modulation patterns do not matter, and we 
may represent -B(x) by a (finite) collection of point sources 

fi(x) =5^fifc(5(x-Xfc) (9) 

k 

X see Appendix 
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with the agreement that Bf^ > and that the spacing between the is not less 
than the finest resolvable scale Iq ~ min(|kj|^^). Since we do not wish to introduce any 
bias into the brightness distribution, except for its localization in the prior region, we 
set dP{B) = Ylk^i'^k) d^k where ^(x) is a pdf concentrated in the prior region. This 
definition of dP{B) is insensitive to the amplitude of B but sensitive to its support. 
Applying £-q reduces now to an elementary calculation, and Equation (jHJ becomes 



r 



CP-'){M,) (10) 



with 



(M^) = jjt I d^i{-K)Mi^^){-^,t) (11) 



(M^M,)= j dt j dt' [d^^{x)M,^^){x,t)Mi^,){x,t') (12) 

n 

7 =l-p-'J2Bl (14) 

n 

Equation ()10j) is our main result. The parameter < 7 < 1 is a 'filling factor': 
7 = corresponds to a single unresolved source, whereas 7 — > 1 stands for a prior 
region densely covered with sources. If the diameter of the prior region exceeds the 
angular pitch, then 7 — > 1 indicates the loss of modulation, while 7 = retains full 
modulation, yet at an unknown phase. Note that (3 is the quantity which actually is to 
be estimated. Its occurrence in Equation ()10|) is, however, unproblematic. On the right 
hand side, (3 is absorbed in the Lagrange multiplier (. On the left hand side, (3 occurs 
in the regularization only, where it may be replaced by a simpler estimate Bq (Sect. 14. ip 
without qualitatively changing the solution. Since ( only affects the scaling of the right 
hand side vector in Equation (fTUjl . its adjustment for unbiasedness (Eq. ^ is equivalent 
to re-scaling w — wr(w ■ (M))^^, which is easily implemented numerically. 

It remains to select the pdf ^(x). Our choice is largely ad hoc, and must be justified 
by demonstrating that the resulting estimator b remains efficient and well-behaved even 
if the prior assumptions are violated. This will be done in Sections 14.41 and 14.51 Here we 
quahtatively motivate our choice. First of all, ^(x) should be simple and depend on not 
more than a few parameters in order to facilitate operational data processing. One set 
of parameters which is provided by the RHESSI data products is the centroid xq of the 
HXR emission derived from time-integrated (^ r) RHESSI imaging. It is natural to 
center ^(x) on xq. In addition to the centroid we would like to specify the rough size / of 
the prior region, but not any of its details, since these are not known on time scales as 
small as r. The choice of / refiects the uncertainty of the true scene, and may be based on 
RHESSI imaging or on independent observations. In particular, we may chose I such as 
to cover a solar 'active region' (Howard 1996) derived from magnetic field observations. 
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In any case, / should not be unrealistically small; simulations (Sect. I4.5p suggest that 
I > 10" is a reasonable, conservative, choice. Owing to the inherent rotation symmetry 
of the rotation modulation observing principle, it is suggestive to choose an isotropic 
form for ^(x), although this is by no means a rigorous request. From a practical point of 
view it is also important that the integrals in Equation (jllll and (fT^ can be performed 
analytically as far as possible to save numerical operations. A choice which fulfills the 
above criteria is 

e(x) = (27r/2)-ie'("-"o)'/2;2 . (15) 

Since the coefficients a„j(t) and ^I/j(t) in Equation (P) are slowly varying with time, 
they may be replaced by discrete-time versions a„y and Equations (fTTl [T^ then 
become 



(M^) = ^a^^e-^'"^'''(-)l' f dt cos {m<j),{t)) (16) 

m=0 "^^f 

(M^M,)= J2 Yl ^^^e-^(l™''^(M)P+l"k,Mp) f f dt'x (17) 

s=±lm,n=0 ^ "^^M "'^"^ 

^ ^-l^smnk,aty^,at') cos (m0^(t) + Sn(f)^{t')) 

with 4>^(t) = kj(^)(t) ■ (xo — P(t)) + "^^^ On time scales r <^ Ts/2, mostly terms 
with s=— 1 contribute to Equation (fT7|) . The time integrals in Equations flTB|) and (fT7|) . 
which depend on the actual spacecraft motion, are evaluated numerically. 



3. 3. Choice of time bins 

RHESSI detects individual photons, and their assignment to time bins is an important 
first step of the demodulation procedure. There are two, potentially conflicting, requests 
on the time bins A^. On the one hand, A^ should resolve (say, by a factor 10) the 
instantaneous modulation period 

rr"(t)~T5|k,(t)x(xo-P(t))|-^ [s] (18) 

of a source at position xq. On the other hand, demodulation is only beneficial if there are 
sufficient counts to observe modulation against Poisson noise. This requires, empirically, 
some 5 counts per time bin. In addition, the time bins should be integer fractions of r 
in order to minimize roundoff error. We therefore adopt the following rule: 

A„ = ^ , , , A* = max f — minr/?°x'^(t)') (19) 

1 + round (A;/r) ' V(q(^))' 10 ter ^ V ^ ' 

where (cj) is the average count rate [ct/s] in subcollimator i. The factors 5 and 10 
in Equation (jl9p are empirical. 
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3.4- Numerical implementation 

Equations (fT6|l and (fTTjl are evaluated by an extended trapezoidal rule with intermediate 
step size adapted to the modulation frequency and amplitude. The inversion of 
the matrix on the left hand side of Equation ()10|) is performed by singular value 
decomposition (Golub and Van Loan 1989) with explicit control of condition number. 



4. Discussion 



4.1. Limiting behaviour 

Let us start our discussion by verifying that the estimator (j3]) remains well-behaved and 
meaningful in extreme cases where the prior is pointlike or flat, and where the observed 
flux tends to zero. The limit of infinite flux is discussed in Sect. 14.21 

At very low count rates (<^ 1 ct/r™°'^), modulation is no longer observable, and one 
may thus expect that b should reduce to a simple average, which is an efficient estimator 
for the pure counting noise problem. Indeed, in the hmit /3 — > 0, Equation (|Tnjl yields 
(after adjusting ( to satisfy Eq. ^ 



(20) 



Equation pn|l will be referred to as 'uniform average' since it involves uniform 
weights (w^=const). The last approximation in Equation ()20|) holds if the modulation 
is fast (r™°'^ <^ r) or weak (|kj|/ ^ 1), so that the term with m=l in Equation can 
be neglected. The uniform average 6o provides a simple guess of b, which is - hopefully 
- improved by the more sophisticated estimator We shall see in Sect. 14.41 - 14.5) that 
this is in fact the case. 

The uniform average (j^U)) is also attained if / = or 7 = 1, both representing 
deterministic limits with completely localized or densely filled prior regions. In either 
case, the term (1 - 7)(M^M^) + 7(M^)(M^) reduces to (M^)(M^), so that Equation 
f|Tn|l has the solution tu^ = const, and b reduces to bo. 

Another situation of interest is / — > 00 where the prior becomes globally flat. For 
/ 00 one can derive from Equations ()16II17|1 that (M^) a^^A^ and {M^M^) — > 
aQ^/S.^aouAi/ + \{,ciifj,/S.^Y5fj_y. Thus Equation (fTn|) can be inverted by the Sherman- 
Morrison formula (Press et al 1998), and one finds (after adjusting Q 

-2 1 ^2 



.... ^ where a2 = l + -(l-7)/5A^^. (21) 

a?" 

From Equation (PT|) the uniform average (PUj) is recovered if ^(1 - 7)/5A^^ < 1. 
Otherwise a correction arises which weakly varies with time (Sect. |2I). Bins with large 
modulation amplitude are given less weight in Equation ()21|). which is reasonable 
since their uncertainty is larger. 
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In summary, we have found that the estimator (0)) has well-defined and meaningful 
analytic limits when the flux tends to zero {/3 —>■ 0), and when the prior is completely 
restrictive (/ — ^ 0) or completely nonrestrictive (/ — oo). 

4-2. Geometrical interpretation 

In the limit of inflnite count rates, the counting noise becomes unimportant and 
Equation (fTIH) admits a purely geometrical interpretation. Recalling the assumption 
i?(x, t) ~ B{x.) for t e T, one has that £c\B{b) = J (ixi?(x) $(x) with = w^M^ipL) 
and M^(x) = dt Mi(^^-^{-K,t). Therefore £c\Bip) is a good estimator for h for arbitrary 
B{x) if $(x) is independent of x. The modulation patterns would then be canceled. 
However, the functions M^(x) are neither orthogonal nor complete, so that $(x) 
cannot be made constant by whatever choice of weights w^. Instead, we may try to 
minimize the fluctuations of $(x) within the spatial prior ^(x), and therefore minimize 
/ ^(x)($(x) —const )^ (ix. This yields {M^Mi,)wu oc (M^). Comparing this expression to 
Equation pOj) we see that Equation ()10|) minimizes the fluctuations of $(x) within the 
spatial prior under the maximum-modulating assumption that there is a single point 
source (7 = 1) observed at inflnite count rate (/9 — >• 00). 

The matrix M = {M^My) is generally ill-conditioned, and the weig hts w = M"^(M) 
therefore oscillate. At flnite count rate, the excursions to large \w^\ amplify the Poisson 
noise, and thereby impair the estimator h (Eq. 0]). An optimum estimator must thus 
limit the oscillations of w to a level which commensurates with the Poisson noise. 
Equation ()10|) provides a possible trade-off. 

Fig. 121 illustrates the resulting function $(x) for 7 = and different /3, using aspect 
data of February 26 2004, sub collimators (3,4,5,6), and r = 0.2s. The time bins are 
(2.8, 4.8, 8.5, 15)ms, and the spatial prior (size / = 40") is centered at the source position 
(245", 340") estimated from time-integrated RHESSI imaging. The black mask in Fig. 
121 indicates the 10% level of the spatial prior (Eq. IT^ . In order to assess the constancy 
of $(x) we consider the standard deviation A$ and mean $ inside the black mask. In a 
noiseless world {{3 — 00), the modulation patterns would admit A$/ ($) ~ 0.0058 (Fig. 
12^). Taking the actual counting noise (/3 ~ 60 = 6200 ct/s) into account, A$/($) ~ 
0.024 is still achievable (Fig. [2b)- For comparison, the uniformly weighted (if^=const) 
case is also shown (Fig. |2l3, A$/($) ~ 0.14), which corresponds to the limit /3 — > 0. 
As can be seen, the function $(x) becomes less efficient in canceling the modulation 
patterns when the count rate (3 decreases. But, at the same time, the sensitivity of h 
to Poisson noise increases. The trade-off chosen by the present method is shown in Fig. 

|2l3. 

4-3. Weighting mechanism 

In a generic situation (77^1, /t^O, 6o^1) one may identify three mechanisms by which 
the weights operate. These are most easily discussed in terms of Fourier modes of the 
modulation. First, if the modulation phase is resolved (|kj|Z <^ I) then modulation 
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can 'actively' be countersteered. Secondly, if the modulation phase is not resolved 
but the instantaneous modulation frequency l/r™°*^ (Eq. ITHj) is known, then counts 
with r°^°<^ it T can be suppressed because they do not allow a 'passive' averaging. 
This suppression works for each sub collimator individually. (Although r™°'^ does not 
explicitly show up in Equation it is implicitly contained in the diagonal blocks 

of correlation function (M^Mj,), see Fig. H] below). Thirdly, the matrix (M^M^) 
couples different subcoUimators, which regulates the relative weighting of different 
subcoUimators. 

The different mechanisms are illustrated in Figures El and |31 showing simulated 
data of subcoUimators (4,5,6,7) with xq = (885", 161"), / = 20", 7 = 0.2, and a single 
intense (~ 3x10^ ct/s/subcollimator) point source located at xq. The chosen geometry 
represents a solar limb source. The simulated counts are divided into disjoint intervals 
of duration r = 0.25s (Fig. IH}, and time bins are equal (A^ = 0.0025s) for better 
comparability. The outweighing of modulation is most clearly seen in the coarsest 
subcollimator (Fig. 01 #7) with l/p^ = 0.16, where p7 = 27r/|k7| = 122" is the period 
of subcollimator ^7. The modulation is also - though less efficiently - outweighed in 
subcollimator #6 {l/pe = 0.28), while the finest subcoUimators #5 (//ps = 0.49) and 
#4 (//p4 = 0.85) mostly operate in an averaging mode, with priods of low r^^^ being 
suppressed. 

The different regimes manifest in different forms of the correlation (M^Mj,) (Fig. 
E}. For l/pi <C 1 the matrix {M^Mi,) approximately factorizes into {M^){My) (Fig. HI 
#7). With increasing l/pi, the autocorrelation of a single subcollimator takes the form 
(M(,,i)M(,,i,)) ~ cos(27ri^)e-(*-*')'/2(-^)' with decay time rf ~ {lh)-\\^, ■ Xo)-i (Eq. 

i 

ITTj) . Here, both rf and r™°*^ depend to first order on (t + t')/2, which results in the 
'chirping' behaviour of Fig. |3J subcollimator ^4. The chirp towards low modulation 
frequencies is associated with a drop of the corresponding weights (Fig. |21 top, time 
interval in dashed lines). 

4.4- Performance 

The quality of the demodulation and its robustness against violation of prior 
assumptions have been explored by Monte Carlo simulations. Figure El shows, as an 
example, the relative errors \h — h\/h oi 10^ simulated time intervals r, together with 
uniform values |6 — and the relative Poisson error N^^^ with A^tot the total number 
of counts in r. Dots illustrate a subsample of the simulation; graphs represent the full 
ensemble averages. The spatial prior is centered at xq = (420", -630") and has size I 
= 40". The simulated brightness distribution is a superposition of 10 random sources 
Asexp{ — ^{t — tsY/Tg — |(x — Xs)^//^} with uniform relative amplitudes As/ J2s^s, 
Gaussian positions x^ (mean xq, variance /^/4), uniform sizes Ig € (1" ... 5"), and uniform 
intrinsic time scales G (0.2s... Is). The assumed filling factor is 7 = 0.2, whereas the 
true value scatters from 0.8 to 0.9. Data gaps are neglected. At low count rates, the 
error is dominated by the Poisson noise, and all three types of error coincide. At higher 
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count rates, the error becomes dominated by modulation, and the uniform average no 
longer improves with improving counting statistics. The estimator b performs better, 
and reaches the Poisson limit for integration times r > 8ms and count rates up to 10^ 
ct/s/sub collimator. The latter is, in fact, the practically relevant range - higher count 
rates do not occur because the RHESSI detectors need ^ 8 x 10~^s to recover after each 
photon impact. One may thus be confident that - at least with regard to the (tolerant) 
error measure £{\b — b\/b) - the estimator b works optimally for practical purposes. The 
degradation of\b — b\/b at highest count rates and largest r is attributed to the increased 
susceptibility to prior assumptions, and possibly also to numerical errors as the involved 
matrices become large and weak-conditioned. (Condition numbers up to 10^ occur at 10^ 
ct/s/subcollimator.) Extended simulations including different brightness distributions, 
data gaps, and 7 values yielded similar results. 

4.5. Influence of the prior 

An important practical feature of Bayesian inverse methods is their robustness against 
violation of the prior assumptions. If a wrong prior is used, then the result will be 
degraded but it should not become (much) worse than if no prior was used at all (i.e., if 
the weights were uniform or the prior was flat, / — >■ 00). We shall now demonstrate that 
this goal is met by our demodulation method. Figure El shows simulations performed at 
fixed count rate of 2 ■ 10^ ct/s/subcollimator and fixed integration time r = 0.12s, but 
with varying size / of the prior region, and with different offsets between the true and 
prior centroids. In each simulation, a prior centroid Xq is chosen at random across the 
solar disc, and the true centroid is placed d arc seconds away in random direction. The 
true source has several gaussian components, which are all concentrated in a narrow 
(3") region around the true centroid. The simulated imaging axis P{t) moves within 
the central (<200") region of the solar disc. While the cases d = (0,10,50)" reflect 
realistic uncertainties of RHESSI observations, 50" being a conservative value, the case 
d = 500" is overly pessimistic and is included here for demonstration purposes. The 
relative errors (vertical axis) are defined as in Figure the 'fiat prior' estimator is given 
by Equation (|2T| . Dots again represent a subsample of the simulation, while curves 
represent averages over the full sample (6 ■ 10^). 

If the prior centroid coincides with the true one {d=0" , black solid line), then the 
demodulation reaches the Poisson limit for I ^ 5", i.e., as long as the prior width 
does not exceed the true width by more than a factor ^2. If the true and assumed 
centroids differ [d > 0), then the demodulation does no longer reach the Poisson limit, 
and degradation depends on the ratio d/l. As can be seen, all error curves collapse to 
the optimum {d = 0) one when / ^ d, indicating that discrepancies between true and 
prior centroids are tolerated up to the size of the prior region. Except for very large and 
unrecognized prior errors {d=500" , I < 20"), the demodulation performs better than 
the uniform average, and approaches the fiat prior estimate as Z ^ 00. The flat prior, 
in turn, performs better than the uniform average. The solar diameter is 1920" , so that 
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/ > 2000" effectively represents a trivial prior, which only indicates that the photons 
came from the sun. 

As a practical implication, we learn from Figure El (and from similar simulations) 
that, for realistic uncertainties (say, d < 50"), demodulation is beneficial compared to 
the uniform average if I ^ 10", and it is beneficial compared to the fiat prior if I > 20". 

4-6. Example of a solar eruption 

A recent solar eruption, to which also Figures ^]21 refer, occurred on February 26, 2004. 
RHESSI has observed the whole eruption, and Figure shows a part of the impulsive 
rise phase, where most temporal fine structures are expected. The observed counts of 
sub collimators (1,3,4,5,6,7,8,9) are divided into disjoint intervals of duration r = 0.12s, 
while subcollimator 2 is rejected due to increased background (Smith et al. 2002). The 
centroid Xo = (245", 340") and size I = 30" of the spatial prior are taken from a long- 
exposure (10s) RHESSI imaging and in agreement with the simulations of Sect. 14.51 
The filling factor is assumed to be 7 = 0.4. The average count rate is about 4000 
ct/s/subcollimator, and energies from 5 to 15 keV are used. 

In order to remove some arbitrariness of the time binning and to assess the 
quality of the demodulation b we proceed as follows. By assumption, the true scene 
is approximated as piecewise constant in time. If this assumption was true then a 
second demodulation b^h with intervals of equal duration r but shifted by r/2 should 
yield a similar result. By comparing b with bsh we may thus gain an estimate on the 
accuracy of the demodulation, and by considering 6avg = |(^ + ^sh) we may remove some 
arbitrariness of the time binning. 

Figure [7| (top panel) shows the estimator 6avg, together with the uniform average 60 
as the simplest possible guess. Panel b) shows the discrepancy between b and its time- 
shifted version 6sh- As can be seen, the relative discrepancy \b — 6sh|/&avg is throughout 



small (6%). The pure Poisson error ^J^^c^w'j^ is also shown for comparison (gray 

line). The residuals b — bsh exceed the pure Poisson error, and this excess is due to 
uncertainties of the weights w. The latter are caused by the principal reasons discussed 
in Sect. 14.21 but also have contributions from the uncertainty of the prior centroid 
Xo and -size /, as well as from the approximate form of Equation and its energy- 
averaged coefficients {ani{t) , ■ipi{t)} , possible errors of the aspect data {P(t), argkj(t)}, 
and from violation of the piecewise constant-in-time approximation of the true scene. A 
complete disentangling of different sources of errors is difficult and somewhat speculative. 
However, we may empirically compare the residuals b — bsh to the residuals 6avg — ^fi 
(Fig. IZh) and &avg — &o (Fig. Eli). This shows a clear order of residual amplitudes b) < c) 
< d), in agreement with the simulations. We may thus be confident that the fiat prior 
estimate improves the uniform average, and that the demodulation improves the fiat 
prior estimate. All residuals are centered about zero, in agreement with unbiasedness. 
By comparing simulation results with the residuals of Fig. [3b-d), and with famihes 
of similar real-data demodulations with varying r (not shown), we conclude that the 
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demodulation error is in the order of the residuals b — bsh (Fig. Ujp). As may be seen 
from Figure [7^), many of the excursions of bo - especially towards low count rates (data 
gaps) - are absent in 6avg- They are therefore most likely to be instrumental. Not 
all of the data gaps can, however, be removed by the demodulation: see, e.g., before 
01:53:50. Here, the collective dropout of several detectors during > r inhibits successful 
compensation. 

5. Summary 

We have developed a unbiased linear Bayes estimator for the photons arriving in front 
of the RHESSI optics, which applies in situations where imaging information is less 
in demand than (spatially integrated) temporal evolution. The prior assumptions 
involve time-independence of the true brightness distribution during r -C T5, and 
an a gaussian a priori pdf for the source density on the solar disc. The estimator 
minimizes the expected quadratic deviation of true and retrieved unmodulated counts, 
while enforcing agreement of their expectation values. Non-overlapping time intervals 
r are independent. Geometrically, the algorithm tries to cancel the spatial transmission 
patterns (modulation patterns) of the RHESSI optics by a suitable linear combination 
of patterns belonging to the time interval r. The degree to which canceling is beneficial 
depends on the counting noise, and the algorithm constructs a trade-off between Poisson- 
and modulational uncertainties of the estimator. Monte Carlo simulations show that 
the mean relative error of the demodulation reaches the Poisson limit, and demonstrate 
robustness against violation of the prior assumptions. An application to a solar eruption 
is also discussed. 

The present method is limited in several ways. First, any non-solar background is 
neglected. Secondly, the use of sharp time intervals r brings along the computational 
advantage that the data can be split and the results merged in the end, but at the cost 
of possible artifacts at interval boundaries. The use of larger and smoothly tapered time 
intervals would help, but is computationally demanding. 

6. Appendix: Condition number inequalities 

The condition number of a positive definite matrix is defined as the ratio of its 
largest eigenvalue to its smallest eigenvalue. We first ask for simple bounds on the 
condition number of the matrix A^,^ = + \^5^u with positive A^. Since A is 

symmetric all its eigenvalues ly between the minimum and maximum of the Rayleigh 
quotient R = (x. Ax) = (x ■ A)^ + A^x^ with |xp = 1 (Euclidian norm). The 
minumum of R is bounded by i? > min(^^A^x^) = min(A^), while the maximum is 
bounded by i? < |xp|Ap + max(^^ A^x^) = |Ap -t- max(A^). Therefore, cond(A) < 
(|Ap -Fmax(A^))/min(A^). 

Next we consider the matrix £^bA. Since it represents an average over different A 
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with different principal directions, we may expect that 

cond(^BA) < maxcond(A) . (22) 

B 

To show that this is true it suffices to show that for any two positive definite 
matrices Ai and A2 of equal size the inequality 

cond(Ai + A2) < max (cond(Ai), cond(A2)) (23) 

holds. Equation (j^^ then follows by induction over i with Ai = dP{Bi)l\[Bi). 
(Assuming that the true scenes may be labeled by discrete labels i. We shall not prove 
this assertion, but call it plausible in view of the finite resolution of the modulation 
patterns.) Equation is then easily verified by direct calculation: 



cond(Ai + A2 



< 



max((x, Aix) + (x, A2X)) 
min((x, Aik) + (x, A2X)) 
max(x, Aix) + max(x, A2X) 
min(x, Aik) + min(x, A2X) 
cond(Ai) min(x, Aix) + cond(A2) min(x, A2X) 
min(x, Aix) + min(x, A2X) 
< max(cond(Ai), cond(A2)) . 

The last line follows by either assuming cond(Ai)<cond(A2) or cond(Ai)>cond(A2). 

The author thanks G. Hurford, A. Csillaghy, R. Schwartz, A. Benz, C. Wigger, P. 
StHilaire and E. Kirk for helpful discussions and comments. 
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01:53:38 01:53:40 01:53:42 01:53:44 01:53:46 
Start Time (26-Feb-04 01:53:37) 



Figure 1. Rotation modulation observing principle. Top: the probability that a 
photon emitted by the source (S) reaches the detector (D) is diminished if its path 
(dotted) penetrates the bars (black) of the grids (Gi,G2). As the grids rotate, 
the source is periodically shadowed and released, so that the observed flux exhibits 
characteristic modulations (bottom) , from which the true scene in the solar plane (x,y) 
can be recovered. The RHESSI instrument has 9 pairs of grids ('subcoUimators'), of 
which only #6 is shown. Times with zero count rates represent data gaps. 
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Figure 2. The function $(x) obtained from subcoUimators (3,4,5,6) and time 
interval t=0.2s starting at Feb. 26 2004, 01:53:39.850 UT. a) for infinite count rate 
(A$/$=0.0058); b) adapted to the actual count rate (A$/$=0.024); c) in the hmit 
of zero count rate (A$/$=0.14). For A$/$ — > the effect of modulation is canceled. 
See text. 
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Figure 4. Correlation {M^M^) for 2.75s < t < 3s in Fig. |21 Diagonal blocks contain 
the temporal autocorrelations of the 4 subcoUimators, with time running from left to 
right. 
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Figure 5. Simulated relative errors of demodulation, uniform average, and the Poisson 
limit. The curves represent averages over 10^ samples, some of which are shown as 
dots. All subcollimators are used, 7 = 0.1, and I = 40". A color rendering of this 
figure is available in the online version. 
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Figure 6. Influence of the prior on tlie demodulation. The true source is d arc seconds 
away from the prior centroid. The simulation explores different prior sizes I at fixed 
count rate 2 • 10* ct/s/subcoUimator and time interval r — 0.12s, and displays relative 
errors similar as in Fig. The 'demodulation' is given by Eq. 1)10(1 : the 'uniform 
average' by Eq. H2()|l : the 'flat prior' by Eq. pi(l . A color rendering of this figure is 
available in the online version. 
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Figure 7. Real-data application of demodulation, a) demodulation (solid) and 
uniform average (dotted) of subcoUimators (1,3,4,5,6,7,8,9) and r = 0.12s. The 
demodulation &avg is an average over two solutions b and bsh (Eqns. ^ 1101) which 
differ only by an offset r/2 of the time intervals, b) discrepancy of the two solutions, 
together with the pure Poisson error (gray line) . c) residuals between &avg and the flat 
prior estimate ba (Eq. I21f) . d) residuals between &avg and the uniform average bo (Eq. 



